home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Business Shareware
/
Business Shareware.iso
/
start
/
hobby
/
aa_51
/
omercury.c
< prev
next >
Wrap
C/C++ Source or Header
|
1993-02-13
|
2KB
|
90 lines
/* Orbital elements and perturbations for the planet Mercury
* using formulas given by Meeus.
* All values computed are for equinox of date.
*/
#include "planet.h"
#include "kep.h"
/* Calculate orbital elements for Mercury
*/
int omercury(e,J)
struct orbit *e;
double J;
{
/* Orbital elements referred to equinox of date: */
e->epoch = J;
/* calculate all the mean anomalies */
manoms(J);
e->M = M1;
/* semimajor axis */
e->a = 0.3870986;
/* eccentricity */
e->ecc = ( -0.000000030*T + 0.00002046)*T + 0.20561421;
#if OFDATE
e->equinox = J;
/* inclination */
e->i = ( -0.0000183*T + 0.0018608)*T + 7.002881;
/* argument of the perihelion */
f = ( 0.0001208*T + 0.3702806)*T + 28.753753;
e->w = mod360(f);
/* longitude of the ascending node */
f = ( 0.0001739*T + 1.1852083)*T + 47.145944;
e->W = mod360(f);
/* mean longitude */
f = ( 0.0003011*T + 149474.07078)*T + 178.179078;
e->L = mod360(f);
#else
e->equinox = J2000;
/* inclination */
e->i = ((-3.5e-8*T + 6.9e-7)*T - 0.0059556)*T + 7.010678;
/* argument of the perihelion */
f = ((4.3e-8*T + 7.445e-5)*T + 0.2842765)*T + 28.839814;
e->w = mod360(f);
/* longitude of the ascending node */
f = ((-6.8e-8*T - 8.844e-5)*T - 0.1254715)*T + 48.456876;
e->W = mod360(f);
/* mean longitude */
f = e->M + e->w + e->W;
e->L = mod360(f);
#endif
return(0);
}
/* corrections due to perturbations
* These are added in after solving Kepler's equation.
*/
int cmercury(e)
struct orbit *e;
{
double p;
/* Mercury's perturbations in longitude: */
f = (5.0*M2 - 2.0*M1 + 12.220)*DTR;
p = 0.00204*cos(f);
f = (2.0*M2 - M1 - 160.692)*DTR;
p += 0.00103*cos(f);
f = (2.0*M5 - M1 - 37.003)*DTR;
p += 0.00091*cos(f);
f = (5.0*M2 - 3*M1 + 10.137)*DTR;
p += 0.00078*cos(f);
e->L += p*DTR;
/* Mercury's perturbations in radius vector */
f = (2.0*M5 - M1 + 53.013)*DTR;
p = 0.000007525*cos(f);
f = (5.0*M2 - 3.0*M1 - 259.918)*DTR;
p += 0.000006802*cos(f);
f = (2.0*M2 - 2.0*M1 - 71.188)*DTR;
p += 0.000005457*cos(f);
f = (5.0*M2 - M1 - 77.75)*DTR;
p += 0.000003569*cos(f);
e->r += p;
return(0);
}